--with profile
--without type_check
--power((((power((27+power(69,1/3))/2,1/3))/3) + (power((2/(27+))),1/3)) , 30000)

--****
-- converted documentation to eudoc markup and some tests by ne1
-- simply distributed with euqt
--
-- euphoria library for handling large fixed integer numbers
--{{{
--+-------------------------+--------------+--------+
--| Lucius Lamar Hilley III | ICQ: 9638898 | Unkmar |
--+-----------------+-------+ AIM: LLHIII  |--------+
--| bigfixedmath.e  |       +--------------+        |
--| Version 0.67.02 | http://www.unkmar.com/        |
--|   2008/06/08    | http://www.hilleyonline.com/  |
--+-----------------+-------------------------------+
--}}}
--
--default namespace bfi  base = 10  places = 40


namespace bfi


--****
-- Repair Road Map
--
-- * [ ] update docs, each step of the way.
-- * [ ] create and maintain bfn_sanity.exw
-- * [ ] create comments through out the file as you figure things out again.
-- * [x] remove log10()
-- * [ ] create bigfixednumber() type
-- * [ ] fix places in cmp, add, sub, mul, div
-- * [ ] create db2b, and b2b.
-- * [x] rename plain_power() as ipwr()
-- * [ ] split integer portion and decimal portion of pwr()
-- * [ ] then use ipwr() and exp()
-- * [x] rename is_textnum() to textnumber()
-- * [ ] make textnumber() more flexible, allow for spaces and commas
-- * [x] fix cmp() to handle negative numbers
-- * [ ] fix cmp() to be places aware


public constant VERSION = 6704
constant LN_LIMIT = 1.1 --used in the natural log() function

-- MAX_INTEGER = power(2, 53) -- hidden bit - http://en.wikipedia.org/wiki/IEEE_754
public constant MAX_BASE = floor(sqrt(power(2, 53))) -- square root of MAX_INTEGER
-- MAX_BASE = power(5 * 683 * 27791, 2)
-- MAX_BASE = 94906265                       --                    94906265
-- MAX_FACTOR = 2*3*(3*3)*5*7*11*13*17*19*23 -- = 9699690 *(3*3) = 87297210

integer base, places, bplaces

object void, temp
void = 1
--~ temp = floor(remainder(VERSION, {1000000, 10000, 100})/{10000, 100, 1})
--~ printf(1, "bigfixedmath.e\n  Version: %d.%02d.%02d\n", temp)

public type flatseq(object x)
  -- check for flatness.
  -- Any embedded sequences?
  if atom(x) then
    return 0
  end if
  for index = 1 to length(x) do
    if sequence(x[index]) then
      return 0
    end if
  end for

  return 1
end type

-- test for text number: "-1,234.567 890 123"
-- returns 1 for true and 0 for false
public type textnumber(object x)
  integer i
  sequence temp

  -- check for flatness.
  -- Any embedded sequences?
  if not flatseq(x)then
    return 0
  end if

  temp =  (x = '.')
  temp += (x = ',')
  temp += (x = ' ')
  temp += (x = '-')
  temp += (x = '+')
  temp += ((x >= '0') and ('9' >= x ))

  i = find(0, temp)
  i += find('.', x[(find('.', x) + 1)..length(x)])
  i += find('-', x[2..length(x)])
  i += find('+', x[2..length(x)])

  return not i
end type

--**
-- [Loose test] - all values in flatseq() should be less than base.
  -- bigfixednumber has the following structure
  -- {base,    number,  places}
  -- {integer, flatseq, integer}
public type bigfixednumber(object x)
  if sequence(x) and (length(x) = 3) and integer(x[1]) and flatseq(x[2]) and integer(x[3]) then
    return 1
  else 
    return 0
  end if
end type

--**
--Removes leading values.
--Example: trim_head("0000.04", '0') = ".04"
public function trim_head(sequence s, atom a)
  integer f

  f = 1
  while (length(s) > f and s[f] = a) do
    f += 1
  end while

  return s[f..length(s)]
end function

--**
--Removes trailing values.
--Example: trim_tail("123.000", '0') = "123"
public function trim_tail(sequence s, atom a)
  integer f

  f = length(s)
  while (f > 1 and s[f] = a) do
    f -= 1
  end while
  if s[f] = '.' then
    f -= 1
  end if

  return s[1..f]
end function

--**
--bigfixedmath number compare function.
--results similar to compare().
public function cmp(bigfixednumber bn1, bigfixednumber bn2)
  atom result
  integer neg, ln1, ln2, b1, b2
  sequence n1, n2
  
  b1 = bn1[1]
  b2 = bn2[1]

  n1 = bn1[2]
  n2 = bn2[2]

  neg = 0
  if (0 > n1[1]) then
    neg -= 1
    n1[1] *= -1
  end if

  if (0 > n2[1]) then
    neg += 1
    n2[1] *= -1
  end if
  
  if (neg) then
    result = neg
  else
    ln1 = length(n1)
    ln2 = length(n2)
    if (ln1 > ln2) then
      n2 = repeat(0, ln1 - ln2) & n2
    elsif (ln2 > ln1) then
      n1 = repeat(0, ln2 - ln1) & n1
    end if
    result = compare(n1, n2)

    if (0 > bn1[2][1]) then
      result = 0 - result
    end if
  end if
  
  return result  
end function

--**
-- unsigned integer addition for big numbers --
public function uadd(flatseq n1, flatseq n2)
  integer ln1, ln2, lb
  sequence result

  ln1 = length(n1)
  ln2 = length(n2)
  if (ln1 = ln2) then
    lb = ln1
    result = n1 + n2
  elsif (ln1 > ln2) then
    lb = ln1
    result = n1
    result[(ln1 - ln2 + 1)..ln1] += n2
  else
    lb = ln2
    result = n2
    result[(ln2 - ln1 + 1)..ln2] += n1
  end if
  
  for A = lb to 2 by -1 do
    while (result[A] >= base) do
      result[A] -= base
      result[A-1] += 1
    end while
  end for

  if (result[1] >= base) then
    result = uadd({0} & result, {0})
  end if

  return result
end function

--**
-- unsigned integer difference for big numbers --
--   generally used for subtraction            --
public function udiff(flatseq n1, flatseq n2)
  integer ln1, ln2, lb
  sequence result

  ln1 = length(n1)
  ln2 = length(n2)
  if (ln1 = ln2) then
    lb = ln1
    if (compare(n1, n2) = 1) then
      result = n1 - n2
    else
      result = n2 - n1
    end if
  elsif (ln1 > ln2) then
    lb = ln1
    result = n1
    result[(ln1 - ln2 + 1)..ln1] -= n2
  else
    lb = ln2
    result = n2
    result[(ln2 - ln1 + 1)..ln2] -= n1
  end if
  
  for A = lb to 2 by -1 do
    while (0 > result[A]) do
      result[A] += base
      result[A-1] -= 1
    end while
  end for
  
  lb = length(result)
  if (lb > 1) then
    for A = 1 to lb do
      ln1 = A
      if (result[A] != 0) then
        exit
      end if
    end for
    result = result[ln1..lb]
  end if

  return result
end function

--**
-- unsigned integer subtraction for big numbers --
public function usub(sequence n1, sequence n2)
  sequence result

  result = udiff(n1, n2)
  if (length(n2) > length(n1)) then
    result[1] *= -1
  elsif (length(n2) = length(n1)) and (compare(n2, n1) = 1) then
    result[1] *= -1
  end if
  
  return result
end function

--**
-- unsigned integer multiplication for big numbers --
public function umul(sequence n1, sequence n2)
--I should slice for multiply instead 
  integer ln1, ln2, lb
  sequence result

  ln1 = length(n1)
  ln2 = length(n2)
  if (ln1 >= ln2) then
    lb = 1
    result = n1 * n2[ln2]
    for A = ln2-1 to 1 by -1 do
      result = uadd(result, (n1 * n2[A]) & repeat(0, lb)) --slice instead of repeat
      lb += 1
    end for
  else
    lb = 1
    result = n2 * n1[ln1]
    for A = ln1-1 to 1 by -1 do
      result = uadd(result, (n2 * n1[A]) & repeat(0, lb)) --slice instead of repeat
      lb += 1
    end for
  end if
  
  for A = length(result) to 2 by -1 do
    while (result[A] >= base) do
      result[A] -= base
      result[A-1] += 1
    end while
  end for

  if (result[1] >= base) then
    result = uadd({0} & result, {0})
  elsif (result[1] = 0) then
    result = usub(result, {0})
  end if

  return result
end function

--**
-- unsigned integer division for big numbers
public function udiv(sequence n1, sequence n2)
  integer c, pos, divisor, quotient, ln1, ln2
  sequence result, dividend, nn1, nn2
  object temp
  
  if equal(n2, {0}) then -- x/0 - attempt divide by zero
    ? 1/0
  end if

  ln1 = length(n1)
  ln2 = length(n2)

  while (ln1 > 1) and (ln2 > 1) do
    if (n1[ln1] != 0) or (n2[ln2] != 0) then
      exit
    end if
    ln1 -= 1
    ln2 -= 1
  end while    

  nn1 = n1[1..ln1]
  nn2 = n2[1..ln2]
  
  c = cmp({base, nn1, places}, {base, nn2, places})
  if (c = 0) then -- x/x - divide by self
    return {1}
  elsif (c = -1) then -- (x - 1)/x - result less than 1.
    return {0}
  end if
  
  pos = length(nn2) - 1
  dividend = nn1[1..pos]

  divisor = nn2[1]
  if (length(nn2) > 1) then
    divisor *= base
    divisor += nn2[2]
  end if

  result = nn1 * 0
  while (length(nn1) > pos) do
    pos += 1
    dividend &= nn1[pos]
    while (cmp({base, nn2, places}, {base, dividend, places}) = 1) and (length(nn1) > pos) do
      pos += 1
      dividend &= nn1[pos]
      dividend = usub(dividend, {0})
    end while
    dividend = usub(dividend, {0})

    if (length(dividend) > 1) then
      temp = dividend[1]
      temp *= base
      temp += dividend[2]
    else
      temp = dividend[1]
    end if

    quotient = floor(temp / divisor)
    temp = umul({quotient}, nn2)
    while (cmp({base, dividend, places}, {base, temp, places}) > -1) do
      quotient += 1
      temp = uadd(temp, nn2)
    end while
    quotient -= 1
    temp = usub(temp, nn2)
    dividend = usub(dividend, temp)
    result[pos] = quotient
  end while

  pos = 1
  while (result[pos] = 0) do
    pos += 1
  end while
  result = uadd(result[pos..length(result)], {0})

  return result
end function

--**
-- Conversion of text numbers to base10
public function b10(sequence s)
  integer f, lr, neg
  sequence before_point, after_point, result

  result = s

  if atom(s[1]) then

    neg = 0
    if (result[1] = '-') then
      result = result[2..length(result)]
      neg = 1
    elsif (result[1] = '+') then
      result = result[2..length(result)]
    end if

    -- remove commas
    lr = length(result)
    f = find(',', result)
    while (f) do
      result[1..lr-1] = result[1..f-1] & result[f+1..lr]
      f = find(',', result)
      lr -= 1
    end while
    result = result[1..lr]

    -- remove spaces
    lr = length(result)
    f = find(' ', result)
    while (f) do
      result[1..lr-1] = result[1..f-1] & result[f+1..lr]
      f = find(' ', result)
      lr -= 1
    end while
    result = result[1..lr]

    f = find('.', result)
    result -= '0'
    before_point = ""
    after_point = ""
    if (f) then
      before_point = result[1..f-1]
      after_point = result[f+1..lr]
    else
      before_point = result
      after_point = repeat(0, places)
    end if

    lr = length(after_point)
    if (lr > places) then
      after_point = after_point[1..places]
    elsif (places > lr) then
      after_point &= repeat(0, places - lr)
    end if
    result = before_point & after_point
    result = trim_head(result, 0)

    if neg then
      result[1] *= -1
    end if
    
    return {10, result, places}
  else
    for A = 1 to length(s) do
      result[A] = b10(s[A])
    end for

    return result
  end if
end function

--**
-- big number Base to Base conversion
--   Converts a big number           
--   to the specified base in tb     
-- ;Generally internal use only       

public function ib2b(sequence bfn, atom to_base)
  integer current_base, neg, num_base --, num_place
  sequence num, num2, result

  if atom(bfn[1]) then
    num_base = bfn[1]
    num  = bfn[2]
--    num_place = bfn[3]
  
    neg = (0 > num[1])
    num2 = num
    if (neg) then
      num2[1] *= -1
    end if
  
    current_base = base
    base = to_base
    result = uadd({num2[1]}, {0})
    for A = 2 to length(num2) do
      result *= num_base
      result = uadd(result, {num2[A]})
    end for
    base = current_base
  
    if (neg) then
      result[1] *= -1
    end if

    result = {to_base, result}
  else
    result = bfn
    for A = 1 to length(bfn) do
      result[A] = ib2b(bfn[A], to_base)
    end for
  end if

  return result
end function

--**
-- integer addition for big numbers
public function iadd(sequence bn1, sequence bn2)
  integer neg, b, b1, b2
  sequence result, n1, n2

  b1 = bn1[1]
  b2 = bn2[1]

  n1 = bn1[2]
  n2 = bn2[2]

  b = base
  base = b1
  if (b1 != b2) then
    n2 = ib2b(bn2, b1)
    n2 = n2[2]
  end if

  neg = 0
  if (0 > n1[1]) then
    neg += 1
    n1[1] *= -1
  end if

  if (0 > n2[1]) then
    neg += 2
    n2[1] *= -1
  end if

  if (neg = 1) then
    result = usub(n2, n1)
  elsif (neg = 2) then
    result = usub(n1, n2)
  else
    result = uadd(n1, n2)
    if (neg) then
      result[1] *= -1
    end if
  end if
  
  base = b
  return {b1, result}
end function

--**
-- integer subtraction for big numbers
public function isub(sequence bn1, sequence bn2)
  sequence temp

  temp = bn2[2]
  temp[1] *= -1

  return iadd(bn1, {bn2[1], temp})
end function

--**
-- integer multiplication for big numbers
public function imul(sequence bn1, sequence bn2)
  integer b, b1, b2, neg
  sequence n1, n2

  b1 = bn1[1]
  b2 = bn2[1]

  n1 = bn1[2]
  n2 = bn2[2]

  if equal(n1, {0}) or equal(n2, {0}) then
    return {b1, {0}}
  end if

  b = base
  base = b1
  if (b1 != b2) then
    n2 = ib2b(bn2, b1)
    n2 = n2[2]
  end if

  neg = 0
  if (0 > n1[1]) then
    neg += 1
    n1[1] *= -1
  end if

  if (0 > n2[1]) then
    neg -= 1
    n2[1] *= -1
  end if

  n1 = umul(n1, n2)
  if (neg) then
    n1[1] *= -1
  end if

  base = b
  return {b1, n1}
end function

--**
-- integer division for big numbers
public function idiv(sequence bn1, sequence bn2)
  integer b, b1, b2, neg
  sequence n1, n2

  b1 = bn1[1]
  b2 = bn2[1]

  n1 = bn1[2]
  n2 = bn2[2]

  if equal(n2, {0}) then
    return 1/0
  end if

  if equal(n1, {0}) then
    return {b1, {0}}
  end if

  if equal(n1, n2) then
    return {b1, {1}}
  end if

  b = base
  base = b1
  if (b1 != b2) then
    n2 = ib2b(bn2, b1)
    n2 = n2[2]
  end if

  neg = 0
  if (0 > n1[1]) then
    neg += 1
    n1[1] *= -1
  end if

  if (0 > n2[1]) then
    neg -= 1
    n2[1] *= -1
  end if

  n1 = udiv(n1, n2)
  if (neg) then
    n1[1] *= -1
  end if

  base = b
  return {b1, n1}
end function

--**
-- convert text number to big number

public function text_num(sequence num)
  sequence result

  result = b10(num)
  if (base != 10) then
    result = ib2b(result, base)
  end if

  return result
end function

--**
-- convert big number to text number
public function num_text(sequence s)
  atom Base
  integer neg, lr
  sequence result, num

  if atom(s[1]) then
    Base = s[1]
    num  = s[2]
    neg = (0 > num[1])
    if (Base != 10) then
      result = ib2b(s, 10)
      result = result[2]
    else
      result = num
    end if

    if (neg) then
      result[1] *= -1
    end if
    result += '0'

    if (places) then
      lr = places - length(result) + 1
      if lr > 0 then
        result = repeat('0', lr) & result
      end if
      lr = length(result)
      result = result[1..lr-places] & "." & result[lr-places+1..lr]
    end if

    if (neg) then
      result = "-" & result
    end if

  else
    result = s
    for A = 1 to length(s) do
      result[A] = num_text(s[A])
    end for
  end if

  return result  
end function

--**
-- bignum add() function
public function add(sequence n1, sequence n2)
  integer text_flag
  sequence bn1, bn2

  bn1 = n1
  bn2 = n2 

  text_flag = 0
  if (textnumber(bn1)) then
    text_flag += 1
    bn1 = text_num(bn1)
  end if

  if (textnumber(bn2)) then
    text_flag += 1
    bn2 = text_num(bn2)
  end if  

  bn1 = iadd(bn1, bn2)

  if (text_flag = 2) then
    bn1 = num_text(bn1)
    if (find('.', bn1)) then
      bn1 = trim_tail(bn1, '0')
    end if
  end if

  return bn1
end function

--**
-- bignum sub() function
public function sub(sequence n1, sequence n2)
  integer text_flag
  sequence bn1, bn2

  bn1 = n1
  bn2 = n2 

  text_flag = 0
  if (textnumber(bn1)) then
    text_flag += 1
    bn1 = text_num(bn1)
  end if

  if (textnumber(bn2)) then
    text_flag += 1
    bn2 = text_num(bn2)
  end if  

  bn1 = isub(bn1, bn2)

  if (text_flag = 2) then
    bn1 = num_text(bn1)
    if (find('.', bn1)) then
      bn1 = trim_tail(bn1, '0')
    end if
  end if

  return bn1
end function

--**
-- bignum mul() function
public function mul(sequence n1, sequence n2)
  integer text_flag
  sequence bn1, bn2

  bn1 = n1
  bn2 = n2 

  text_flag = 0
  if (textnumber(bn1)) then
    text_flag += 1
    bn1 = text_num(bn1)
  end if

  if (textnumber(bn2)) then
    text_flag += 1
    bn2 = text_num(bn2)
  end if  

  bn1 = imul(bn1, bn2)
  if (bplaces >= length(bn1[2])) then
    bn1[2] = {0}
  else
    bn1[2] = bn1[2][1..length(bn1[2])-bplaces]
  end if

  if (text_flag = 2) then
    bn1 = num_text(bn1)
    if (find('.', bn1)) then
      bn1 = trim_tail(bn1, '0')
    end if
  end if

  return bn1
end function

--**
-- bignum div() function
public function div(sequence n1, sequence n2)
  integer text_flag
  sequence bn1, bn2

  bn1 = n1
  bn2 = n2 

  text_flag = 0
  if (textnumber(bn1)) then
    text_flag += 1
    bn1 = text_num(bn1)
  end if

  if (textnumber(bn2)) then
    text_flag += 1
    bn2 = text_num(bn2)
  end if  

  bn1[2] &= repeat(0, bplaces)

  bn1 = idiv(bn1, bn2)

  if (text_flag = 2) then
    bn1 = num_text(bn1)
    if (find('.', bn1)) then
      bn1 = trim_tail(bn1, '0')
    end if
  end if

  return bn1
end function

--**
-- routine to set fixed point places
-- set_places(2) is good for dollars
public procedure set_places(integer i)
  integer wanted_places

  wanted_places = i
  if (i) then
    bplaces = floor(log(power(10, wanted_places)) / log(base)) + 1
    places =  floor((log(power(base, bplaces)) / log(10)) + 0.01)
--    ? {base, wplaces, bplaces, places, power(base, bplaces)}
  else
    bplaces = 0
    places = 0
  end if
end procedure

--**
-- set_base() only works for integers
-- It will be fixed.  And soon.
public procedure set_base(integer i)
  base = i
  set_places(places)
end procedure

--This routine is fairly stupid as to value types.
--Is pbase a bignum or a text_num?
--If text_num, then effeciency is loss due to repeated
--Conversion to and from bignum format.
public function ipwr(object pbase, object exponent)
  object result
  atom adj_exponent
  integer text_flag
  sequence pbase_bn, exponent_bn, delta_bn
  sequence Z0, Z1, Z2
  object index, index2

  text_flag = 0
  Z1 = text_num("1")
  if atom(exponent) then
    result = pbase
    adj_exponent = exponent
    if atom(pbase) then
      index = 2
      while (adj_exponent >= index) do
        result *= pbase
        index += 1
      end while
    else
      if textnumber(pbase) then
        pbase_bn = text_num(pbase)
        text_flag = 2
      else
        pbase_bn = pbase
      end if

      result = Z1
      while (adj_exponent >  0) do
        delta_bn = pbase_bn
        index = 1
        index2 = 2
        while (adj_exponent >= index2) do
          index = index2
          delta_bn = mul(delta_bn, delta_bn)
          index2 *= 2 
        end while
        result = mul(result, delta_bn)
        adj_exponent -= index
      end while
    end if
  else
    Z0 = text_num("0")
    Z2 = text_num("2")
    
    pbase_bn = pbase
    exponent_bn = exponent
    if textnumber(pbase) then
      pbase_bn = text_num(pbase)
      text_flag += 1
    end if
    if textnumber(exponent) then
      exponent_bn = text_num(exponent)
      text_flag += 1
    end if

    result = Z1
    while (cmp(exponent_bn, Z0) = 1) do -- (adj_exponent > z0)
      delta_bn = pbase_bn
      index = Z1
      index2 = Z2
      while (cmp(exponent_bn, index2) > -1) do -- (adj_exponent >= index2)
        index = index2
        delta_bn = mul(delta_bn, delta_bn)
        index2 = mul(index2, Z2)
      end while
      result = mul(result, delta_bn)
      exponent_bn = sub(exponent_bn, index)
    end while
  end if
  
  if (text_flag = 2) then
    result = num_text(result)
    if (find('.', result)) then
      result = trim_tail(result, '0')
    end if
  end if
  
  return result
end function

--**
-- Used to get Euler's number
-- e = Euler's number
-- exp() is a power function.
-- e = exp(1)
-- exp(5) = power(e, 5)
public function exp(object x)
  integer neg, text_flag
  object temp, result, old, inc, xi, x1, tt, Z0, Z1
  text_flag = 0

  if atom(x) then
    if (x = 0) then
      return 1
    elsif (0 > x) then
      x1 = -x
      neg = 1
    else
      neg = 0
      x1 = x
    end if
  
    old = 0
    result = 1
    temp = 1
    inc = 1
    xi = x1
    while (old != result) do -- Until it stops changing.
      old = result
      temp *= inc
      result += xi / temp
      xi *= x1
      inc += 1
    end while
  
    if neg then
      return 1/result
    else
      return result
    end if
  else -------------------------------------------------
    if textnumber(x) then
      text_flag = 1
      x = text_num(x)
    end if

    Z0 = text_num("0")
    Z0[2] = {0}
    Z1 = text_num("1")
  
    if (cmp(x, {x[1], {0}}) = 0) then
      return "1"
    elsif (cmp({x[1], {0}}, x) > 0) then
      x1 = sub(Z0, x)
      neg = 1
    else
      neg = 0
      x1 = x
    end if
  
    old = Z0
    result = Z1
    temp = Z1
    inc = Z1
    xi = x1
    while compare(old, result) do -- Until it stops changing.
      old = result
      temp = mul(temp, inc)
      tt = div(xi, temp)
      result = add(result, tt)
      xi = mul(xi, x1)
      inc = add(inc, Z1)
    end while
  
    if (neg) then
      result = div(Z1, result)
    end if

    if (text_flag) then
      result = num_text(result)
      if (find('.', result)) then
        result = trim_tail(result, '0')
      end if
    end if
  end if

  return result
end function

--**
-- This is log using base e
-- Also known as natural log
--public function loge(atom a)
--public function nlog(atom a)
--public function n_log(atom a)
--public function natlog(atom a)
--public function nat_log(atom a)
--public function natural_log(atom a)
public function ln(object a)--(atom a)
--  atom result, current_power, old, a1, ln_limit
--  integer inc
--
--  sequence ln_LIMIT, Z0, Z1
--  sequence result, current_power, old, a1, ln_limit, inc
  object ln_LIMIT, Z0, Z1
  object result, current_power, old, a1, ln_limit, inc
  integer text_flag text_flag = 0

  if atom(a) then
    if 0 > a then
      return -1
    elsif 1 > a then
      return -ln(1/a)
    end if

    if (a > LN_LIMIT) then
      a1 = a
      ln_limit = ln(LN_LIMIT)
      result = 0
      while (a1 > LN_LIMIT) do
        a1 = a1 / LN_LIMIT
        result += ln_limit
      end while
      return result + ln(a1)
    end if
    
  -- x - x^2/2 + x^3/3 - x^4/4 + x^5/5
    a1 = a - 1
    current_power = a1
    inc = 1
    result = current_power
    old = a
    while (result != old) do -- Until it stops changing.
      old = result
      inc += 1
      current_power *= a1
      result -= (current_power/inc)
      inc += 1
      current_power *= a1
      result += (current_power/inc)
    end while
  else ------------------------------------------
    if textnumber(a) then
      text_flag = 1
      a = text_num(a)
    end if

    Z0 = text_num("0")
    Z1 = text_num("1")

    ln_LIMIT = text_num(sprintf("%f", {LN_LIMIT}))
  
    if (cmp(Z0, a) > 0) then
      result = text_num("-1")
    elsif (cmp(Z1, a) > 0) then
      result = sub(Z0, ln(div(Z1, a)))
    elsif (cmp(a, ln_LIMIT) > 0) then
      a1 = a
      ln_limit = ln(ln_LIMIT)
      result = Z0
      while (cmp(a1, ln_LIMIT) > 0) do
        a1 = div(a1, ln_LIMIT)
        result = add(result, ln_limit)
      end while
      result = add(result, ln(a1))
    else
    -- x - x^2/2 + x^3/3 - x^4/4 + x^5/5
      a1 = sub(a, Z1)
      current_power = a1
      inc = Z1
      result = current_power
      old = a
      while compare(result, old) do -- Until it stops changing.
        old = result
        inc = add(inc, Z1)
        current_power = mul(current_power, a1)
        result = sub(result, div(current_power,inc))
        inc = add(inc, Z1)
        current_power = mul(current_power, a1)
        result = add(result, div(current_power,inc))
      end while
    end if

    if (text_flag) then
      result = num_text(result)
      if (find('.', result)) then
        result = trim_tail(result, '0')
      end if
    end if
  end if

  return result
end function

--**
-- bignum's version of the power() function
-- Uses exp() and log() to get its values
-- May return inaccurate results for integer
-- exponents
-- Example: pwr("8", "14")
public function pwr(object pbase, object exponent)
  sequence Z0, result
  integer text_flag text_flag = 0
  -- Derived from R.D. Flowers knowledge
  -- Helpful co-worker

--  if (floor(exponent) = exponent) then -- use ipwr
--    return ipwr(base, exponent)
--  end if
  if atom(pbase) and atom(exponent) then
    if (exponent > 0) then
      return exp(ln(pbase) * exponent)
    else
      return exp(ln(pbase) / (-exponent))
    end if
  else
    if atom(pbase) then
       pbase = text_num(sprintf("%f", {pbase}))
    elsif textnumber(pbase) then
      text_flag = 1
      pbase = text_num(pbase)
    end if
    if atom(exponent) then
      exponent = text_num(sprintf("%f", {exponent}))
    elsif textnumber(exponent) then
      text_flag = 1
      exponent = text_num(exponent)
    end if

    Z0 = text_num("0")
    if (cmp(exponent, Z0) > 0) then
      result = exp(mul(ln(pbase), exponent))
    else
      result = exp(div(ln(pbase), sub(Z0, exponent)))
    end if

    if (text_flag) then
      result = num_text(result)
      if (find('.', result)) then
        result = trim_tail(result, '0')
      end if
    end if
  end if

  return result
end function

--**
-- function for finding the log in any base
-- Examples:
-- where 256 is the power and 2 is the base
-- logb(256, 2) = 8
-- logb(1000, 10) = 3
public function logb(object product, object base)
  sequence result
  integer text_flag text_flag = 0
  if atom(product) and atom(base) then
    return ln(product) / ln(base)
  else
    if atom(product) then
      product = text_num(sprintf("%f", {product}))
    elsif textnumber(product) then
      text_flag = 1
      product = text_num(product)
    end if
    if atom(base) then
      base = text_num(sprintf("%f", {base}))
    elsif textnumber(base) then
      text_flag = 1
      base = text_num(base)
    end if

    result = div(ln(product), ln(base))

    if (text_flag) then
      result = num_text(result)
      if (find('.', result)) then
        result = trim_tail(result, '0')
      end if
      return result
    else
      return result
    end if
  end if
end function

public function root(object product, object exponent)
  sequence result
  integer text_flag text_flag = 0
  if atom(product) and atom(exponent) then
    if 0 > exponent then
      return -1
    else
      return pwr(product, 1/exponent)
    end if
  else
    if atom(product) then
      product = text_num(sprintf("%f", {product}))
    elsif textnumber(product) then
      text_flag = 1
      product = text_num(product)
    end if
    if atom(exponent) then
      exponent = text_num(sprintf("%f", {exponent}))
    elsif textnumber(exponent) then
      text_flag = 1
      exponent = text_num(exponent)
    end if

    if (cmp(text_num("0"), exponent) > 0) then
      result = text_num("-1")
    else
      result = pwr(product, div(text_num("1"), exponent))
    end if

    if text_flag then
      return num_text(result)
    else
      return result
    end if
  end if
end function

-- Don't blame me.  This one was kat's idea. :)
-- 2008-06-04
public function bfn_to_int(sequence bnp, atom bit_size)
  sequence result, bnp2
  integer backward, byte_size

  bnp2 = bnp
  if (bnp2[1] != 256) then
    bnp2 = ib2b(bnp2, 256)
  end if

  byte_size = bit_size/8
  result = repeat(0, byte_size)
  
  backward = length(bnp2[2])
  for forward = 1 to byte_size do
    result[forward] = bnp2[2][backward]
    backward -= 1
    if (backward = 0) then
      exit
    end if
  end for
  
  return result
end function


base = 10
places = 40
bplaces = places

--  set_places(18)
--  temp = text_num({"-1.05", "-1.05", "2"})
--  printf(1, "{%s,\n %s,\n %s\n}\n", num_text(temp))
--  
--  void = add(temp[1], temp[2])
--  --? ib2b(void, 10)
--  void = trim_tail(num_text(void), '0')
--  printf(1, "%s\n", {void})
--  
--  void = sub(temp[1], temp[2])
--  --? ib2b(void, 10)
--  void = trim_tail(num_text(void), '0')
--  
--  printf(1, "%s\n", {void})
--  void = mul(temp[1], temp[2])
--  --void = mul(temp[1], add(temp[3], {11, {2}}))
--  --? ib2b(void, 10)
--  void = trim_tail(num_text(void), '0')
--  
--  printf(1, "%s\n", {void})
--  void = div(temp[1], temp[2])
--  --? ib2b(void, 10)
--  void = trim_tail(num_text(void), '0')
--  
--  printf(1, "%s\n", {void})
--  
--  puts(1, "C:\\>")
--  ? machine_func(26, 0)

-------------------------
--  Dependency Section --
-------------------------
--  x = trim_head(string, char)
--  
--  x = trim_tail(string, char)
--  
--  x = cmp(bn1, bn2)
--  
--  x = uadd(n1, n2)
--    uadd()
--  
--  x = udiff(n1, n2)
--  
--  x = sub(n1, n2)
--    udiff()
--  
--  x = umul(n1, n2)
--    uadd()--carry?
--    usub()--trimming?
--  
--  x = udiv(n1, n2)
--    cmp()
--    usub()--trimming?
--    umul()
--    uadd()
--    usub()
--  
--  x = b10(textnumber)
--    trim_head()
--    b10()
--  
--  x = ib2b(bn, to_base)
--    uadd()--carry?
--    ib2b()
--  
--  x = iadd(bn1, bn2)
--    usub()
--    uadd()
--  
--  x = isub(bn1, bn2)
--    iadd()
--  
--  x = imul(bn1, bn2)
--    ib2b()
--    umul()
--  
--  x = idiv(bn1, bn2)
--    ib2b()
--    udiv()
--  
--  x = add(bn1, bn2)
--    iadd()
--  
--  x = sub(bn1, bn2)
--  
--  x = mul(bn1, bn2)
--    imul()
--  
--  x = div(bn1, bn2)
--    idiv()
--  
--  x = text_num(textnumber)
--    b10()
--    ib2b()
--  
--  x = num_text(bn1)
--    ib2b()
--    num_text()
--  
--  set_places(integer)
--    ln()
--    pwr()
--  
--  set_base(integer)
--    set_places()
--  
--  x = exp(bn1_or_atom)
--    text_num()
--    cmp()
--    sub()
--    mul()
--    div()
--    add()
--  
--  x = ln(bn1_or_atom)
--    text_num()
--    cmp()
--    sub()
--    ln()
--    div()
--    add()
--    mul()
--    
--  x = pwr(bn1_or_atom base, bn1_or_atom exponent)
--    exp()
--    ln()
--    text_num()
--    cmp()
--    mul()
--    div()
--  
--  x = logb(bn1_or_atom product, bn1_or_atom base)
--    ln()
--    div()
--    
--  x = root(bn1_or_atom product, bn1_or_atom exponent)
--    pwr()
--    div()
--    text_num()
--    
--  x = root10(bn1_or_atom product)
--    root()
--    text_num()

